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Abstract 

A new diffuse interface model for a two-phase flow of two incompressible 
fluids with different densities is introduced using methods from rational con- 
tinuum mechanics. The model fulfills local and global dissipation inequalities 
and is also generalized to situations with a soluble species. Using the method 
of matched asymptotic expansions we derive various sharp interface models 
in the limit when the interfacial thickness tends to zero. Depending on the 
scaling of the mobility in the diffusion equation we either derive classical sharp 
interface models or models where bulk or surface diffusion is possible in the 
limit. In the two latter cases the classical Gibbs-Thomson equation has to be 
modified to include kinetic terms. Finally, we show that all sharp interface 
models fulfill natural energy inequalities. 
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1 Introduction 

In recent years diffuse interface models have been successfully used to describe the 
flow of two or more immiscible fluids both for theoretical studies and numerical simu- 
lations. One fundamental advantage of these models is that they are able to describe 
topological transitions like droplet coalescence or droplet break-up in a natural way. 
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In the case of two incompressible, viscous Newtonian fluids the basic diffuse interface 
model is the so-called "Model H", cf. Hohenberg and Halperin [2T] . It leads to the 
Navier-Stokes/ Cahn-Hilliard system 

pd t v + p(v ■ V)v - dw(2r](c)Dv) + Vp = -e div(Vc ® Vc), (1.1) 

divv = 0, (1.2) 

d t c + v • Vc = div(mV/i), (1.3) 

H = e"V(c) - ^Ac. (1.4) 

Here p is the density, v is the mean velocity, Dv = |(Vv + Vv T ), p is the pressure, 
and c is an order parameter related to the concentration of the fluids (e.g. the 
concentration difference or the concentration of one component). Moreover, 77(c) > 
is the viscosity of the mixture, e > is a (small) parameter, which is related to the 
"thickness" of the interfacial region, ip is a homogeneous free energy density and p 
is the chemical potential. Capillary forces due to surface tension are modeled by an 
extra contribution eVcCg) Vc := £:Vc(Vc) T in the stress tensor leading to the term on 
the right-hand side of (11.11) . Moreover, we note that in the modeling diffusion of the 
fluid components is taken into account. Therefore div(mVp) is appearing in (II. 3p . 
where m = m(c) > is the mobility coefficient. 

One of the fundamental modelling assumptions is that the densities of both com- 
ponents as well as the density of the mixture p are constant. Of course, this restricts 
the applicability of the model to situations when density differences are negligible. 
Gurtin et al. [2D] derived this model in the framework of rational continuum mechan- 
ics and showed that it satisfies the second law of thermodynamics in a mechanical 
version based on a local dissipation inequality. 

Lowengrub and Truskinovsky [25] derived a thermodynamically consistent exten- 
sion of the Model H for the case of different densities, which leads to the system: 

pd t \ + p(v • V)v - div S(c, D\) + Vp = -e div(pVc <g> Vc), (1.5) 

d t p + div(pv) = 0, (1.6) 
pd t c + pv ■ Vc = div(ra(c)Vp), (1.7) 

/* = ~P' 2 % + *"V(c) - ~ div(pVc), (1.8) 
oc p 

where S(c,fv) = 2r](c)Dv + A(c)divvl and A(c) is the bulk viscosity coefficient. 
Here the free energy has the density p(e^ 1 ip(c) +e^$-) per unit volume. A simplified 
version of this model has been successfully used for numerical studies, cf. Lee et 
al. [221 S3]- In contrast, there are - to the best of the authors' knowledge - no 
discrete schemes available which are based on the full model (|1.5p -( TL~8"]) . This may 
be due to fundamental new difficulties compared with Model H Ql.ip - fll.4p . For 
instance, the velocity field v is no longer divergence free and the pressure p enters 
the equation for the chemical potential fll.8p . At least analytically, these difficulties 
could be overcome, see Abels [1] for existence of weak solutions. Mathematically 



2 



the coupling of the Navier-Stokes (ll.5p - (ll.6p and the Cahn-Hilliard part (11 . 7j) -f TTTH]) 
is much stronger and the linearized system is very different from the linearization of 
Model H, cf. Abels [2], where strong solutions locally in time are constructed. 

Alternative generalizations of the Model H for the case of different densities were 
presented and discussed by Boyer [5] and Ding et al. [9j. The model by Ding et al. 
consists of fll.ip - fll.4|) . but now for a variable density p = p(c). In order to justify 
this generalization they start from the mass balance equation 

8 tPj + div(p jVj ) = (1.9) 

for the individual fluids j = 1,2 and define the mean velocity v of the mixture as 
volume averaged velocity v = «iVi + u 2 v 2 , where Uj is the volume fraction of fluid j. 
Then flTJ} yields 

divv = 0, (1.10) 

cf. Section [2]below. In contrast to that Lowengrub and Truskinovsky define the mean 
velocity v as mass averaged/bary centric velocity pv = piVi + p 2 v 2 , which yields 

d t p + div(pv) = 0. 

The incompressibility relation (11.101) of course has advantages with respect to nu- 
merical simulations - see the computations related to the model by Ding et al. in 

i- 

Unfortunately, (ll.ip - (ll.4l) seems to be not consistent with thermodynamics in the 
case when p is not constant. Neither global nor local energy inequalities are known 
to hold for (ll.ip - (ll.4l) in that case. The model by Boyer is more complicated. But 
it is derived using a volume averaged mean velocity, which leads to a divergence 
free mean velocity field too. The further derivation of Boyer differs from the one in 
[9] and ours since the starting point are the equations for the conservation of linear 
momentum of each single fluid. Also for this model neither global nor local energy 
inequalities seem to be known, cf. also [I]. 

It is the purpose of the present paper to derive a thermodynamically consistent 
generalization of (ll.ll) - (ll.4p to the case of non- matched densities based on a solenoidal 
velocity field v. More precisely, we will derive the system 

d t (pv) + div(pv <g> v) — div(2^((p)Dv) + Vp = —as div(V<p ® V<p), 

div v = 0, 
d t (p + v ■ V<p = div(m(<p)Vp^), 

together with 



-i / x dp l v l 2 

H v = ae ^ ((f) - aeA(p 



2 ' 

where the order parameter tp = <p 2 — <Pi stands for the difference of volume fractions 
(fj, j = 1, 2^- Here the free energy of the system has the density ae^f^) + ae ^^ 



1 Other choices of order parameters are possible as well - see Section [5] 
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(per unit volume). In comparison with the system derived in [5] there is the additional 
term m ^ ne ec L ua tion for the chemical potential. This term vanishes in the 

case of matched densities, i.e., p = const.. But this term is crucial in the case of non- 
matched densities for consistency with thermodynamics. We note that in contrast to 
the model by Lowengrub and Truskinovsky (ll.5p - fll.8B the usual continuity equation 
(jl.fip is not part of our system. Nevertheless there is conservation of mass in our 
system. More precisely, we have 

dtp + div (^ov - ^y^m(^)V/iJ = 0, (1.11) 

and in Section [2] it will be shown that in fact individual masses are conserved. Note 
that pj is the specific density of fluid j = 1,2 and that p = + P2 ~^ (f. We 

emphasize that according to equation (II. lip the (volume averaged) velocity v does 
not describe the flux of the density. In our model the flux of the density consists of 
the two parts: pv, describing the transport by the mean velocity, and a relative flux 
— P2 ~ pi m((/g)V/x related to diffusion of the components. Hence the diffusion of the 
components relative to the mean velocity leads to a diffusion of the mass density in 
the case that p\ ^ p2- Moreover, we note that in the classical Model H effects related 
to diffusion of the components can play an important role and can lead to Ostwald 
ripening effects and disappearance of small droplets, cf. e.g. pTj . 

The structure of the paper is as follows: In Section [2] we will derive the generaliza- 
tion of Model H, described above in the framework of rational continuum mechanics. 
First we will use a local dissipation inequality and a choice of the energy flux as in 
[2"U] to derive restrictions for the form of the stress tensor and the chemical potential, 
which finally leads to our model after suitable constitutive assumptions. Then we 
briefly discuss the changes in the derivation if the energy flux is not specified at the 
beginning and Liu's Lagrange multiplier method is used, cf. [21]. In Section [3], we 
present a third approach to derive a thermodynamically consistent model, this time 
based on Onsager's variational principle. We consider the more general situation 
when either the system is subjected to gravitational forces or when one additional 
soluble species is present in both fluids. In the latter case, transport effects across 
the interface are taken into account, too, and we derive a diffuse interface analogue 
of Henry's law, cf. Section [3J 

In Section H] we discuss the sharp interface asymptotics in the limit e — > for 
the diffuse interface model together with a soluble species. This is done by using the 
method of formally matched asymptotics. We show that the limit system depends 
essentially on the choice and the scaling of the mobility. Actually, we consider four 
cases related to choosing the mobility degenerate or non-degenerate and letting the 
mobility tend to zero or not. If the mobility m((p) vanishes as e — > 0, we end up with 
the classical model for a two-phase flow with the Young-Laplace law 

-[2r}Dv]±v + \p\tu = (jkv atr(t), 
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where T(t) is the interface between the fluids, u is a unit normal to T(t), k is its mean 
curvature, a is a surface tension coefficient, and [.]! denotes the jump of a quantity 
at F(t) in the direction of v. Moreover, the interface is transported by the velocity 
of the fluid, i.e., 

V-i/-v = atr(t), 

where V is the normal velocity of T(t). If a soluble species with density w is present, 
we obtain the classical Henry condition for the jump of the concentrations of the 
soluble species. 

In the case of a constant mobility, we obtain in the limit e — > 

2(v-i/-V) = mo[Vfi]t-u atr(t), 

for the evolution of the interface and 

[p]±v(y ■ v - V) - [2riDv]±v + \p]±v = okv at T(t), 

for the jump of the stress tensor. Here /i satisfies 

fi = on - ||v| 2 [p]+ - [w]t atr(t), 

fj, is harmonic in the bulk, and m > is a diffusion coefficient related to m. In 
particular the interface is no longer material and diffusion of mass through the bulk 
is still present in the model. In the case of a non-vanishing, degenerate mobility the 
evolution of the interface is governed by the surface diffusion law 

2( v • v - V) = mA r(t) |i, 

where Am) is the Laplace-Beltrami operator of T(t) and m > is a diffusion co- 
efficient related to m. In Section [5] we prove that energy estimates are valid for 
sufficiently smooth solutions of the sharp interface models. Finally, several impor- 
tant identities for the formally matched asymptotics calculations are shown in the 
appendix. 

Acknowledgements: The authors are grateful to Hans Wilhelm Alt for several 
stimulating discussions, which helped to improve the presentation of the derivation 
of the model. This work was supported by the SPP 1506 "Transport Processes at 
Fluidic Interfaces" of German Science Foundation (DFG) through the grants GA 
695/6-1 and GR 1693/5-1. 

2 Derivation of the Model 

In order to derive the diffuse interface model, one assumes a partial mixing of the 
macroscopically immiscible fluids in a thin interfacial region. Therefore let us intro- 
duce some terminology related to mixtures. In the following the fluids are labeled 
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by j = 1,2 and they fill a domain Q C R d . The total mass density of the mixture is 
denoted by p. Moreover, pj denotes the mass density of the fluid j, i.e., 

Mj = / Pj(x) dx 
Jv 

is the mass of the fluid j contained in a set V C f2 and we obtain p = p\ + p 2 . 
Moreover, we denote by Cj = — the mass concentration and note that C\ + c 2 — 1. 

Denoting by Jj the mass flux of fluid j, the mass balance equation in local form 
is given by 

d t pj + div Jj = 0. 

Defining the velocities Vj, j = 1, 2, of the single fluids as Vj = Jj/ pj the mass balance 
equation can be rewritten as 

d t pj + div(pjVj) = 0. 

In what follows we assume that the volume occupied by a given amount of mass of 
the single fluids does not change after mixing, i.e., the excess volume due to mixing 
is zero. If pj is the specific (constant) density of the unmixed fluid j, we introduce 
Uj = j-. The assumption that the excess volume is zero results in 

ui + u 2 =l. (2.1) 
Expressed in terms of the mass concentrations c\ and C2, condition (12. ip reads as 
ci p c 2 p = l 

Pi P2 P Pi P2 

Introducing the mass concentration difference c = c 2 — cj, the above relation implies 
that p = p(c) with a function p defined via 

1 |(c+l) |(1 -c) 



P(c) P2 Pi 

We remark that possible choices for the order parameter in the phase field model 
are the mass concentration difference c, the density difference p := p2 — p\ or the 
difference of volume fractions u := u 2 — u%. 

We now introduce a suitable averaged velocity of the mixture. In contrast to the 
mass averaged/barycentric velocity v given by pv = piVi + P2V2, cf. Lowengrub and 
Truskinovsky (23], we choose the volume averaged velocity v of the mixture as in 
Boyer [5] and Ding et al. jS]. More precisely, we define 

Pi . P2 
V = ItxVi + M 2 V 2 = — Vi + — V 2 , 

Pi P2 
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cf. [51 E]- As a consequence, we obtain, using the fact that the p/s are constants, 
div v = div(| Vl ) + div(gv 2 ) = div(| + f ) = -dt (l + f) = -dtl = 0. (2.2) 

Furthermore, we denote by Jj = Jj — pjV the mass flux of the fluid j relative to the 
velocity v, i.e., 

dtPj + div(pjv) + div Jj = 0. 
Because of (12. 2p and pj = pcj = PjUj, we have 

d t (pcj) + v • V(pcj) + div Jj = 

or equivalently, 

d t Uj + v • Vuj + div Jj = 0, 

where Jj = j-. Because of U\ + u 2 = 1, we require Ji + J 2 = ^ + || = 0, cf. (12. 2p . 
In particular, we obtain 

9 t (pc) + v • V(pc) + div J = 0, (2.3) 

where J = J 2 — Ji- In addition we have 

d t p = d t (pi + p 2 ) — — div(Ji + piv + J 2 + p 2 v) = - div(pv + J x + J 2 ). 

If pi 7^ p 2 , we have in general div(Ji + J 2 ) 7^ 0. Hence the classical continuity 
equation does not hold with respect to the velocity v. This reflects the fact that we 
allow for mass diffusion in the system. 

Instead of c, one could use u := u 2 — ui as order-parameter. Because of 

p = pci + pc 2 = p\Ui + p 2 w 2 , (2.4) 
pc = pc 2 - pci = p 2 u 2 - pi«i (2.5) 

and m 2 = ^y^, u\ = we can also assume that p = p{u) and pc = pc{u). In order 
to have flexibility in the choice of the order parameter, we assume that p = p((p) and 
pc = pc((p) for some order parameter ip such that 

d t pc((p) + v ■ Vpc(^) + div J = 0, (2.6) 

where pc = pc(<^) is the density difference. This is equivalent to 

d(pc) 



dp 



(d t cp + v- V(p) = -div J. (2.7) 



As in Gurtin et al. [20] , we assume that the inertia and kinetic energy due to the 
motion of the fluid relative to the gross motion is negligible. Therefore we consider 
the mixture as a single fluid with velocity v, which satisfies the laws of conservation 
of linear and angular momentum of continuum mechanics with respect to the volume 
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averaged velocity. The density and the stress tensor are assumed to depend on 
additional internal variables like ip and Vy?. I.e., we assume that 

dt(pv) + div(pv <g> v) = divT (2.8) 

for a stress tensor T, which has to be specified by constitutive assumptions. Here 
external forces are neglected for simplicity. 

Because of the constraint div v = 0, the stress tensor T is only determined up to 
pi for some scalar quantity p. Hence we make the ansatz 

T = S-pI where S = S(L>v, y?, Vy?). 

In order to specify the behavior of the mixture, we will make several constitutive 
assumptions on the form of the stress tensor T and the relative mass flux J in the 
following. 

Finally, we assume the relative motion of the fluids to be diffusive, and we intro- 
duce a Helmholtz free energy density f(ip, Vy?) (per unit volume). It will play the 
role of an interfacial energy for the diffuse interface. The total energy in a volume V 
is then obtained as the sum of the kinetic and the free energy, i.e., 

E v (<p,v) = J p(p)^-dx + J f(ip, Vp)dx = J e(v,y9,Vy?) dx, 
where e(v, ip, Vy?) := p(<p)^£- + /(</?, Vy?). 



2.1 Derivation based on a Local Dissipation Inequality and 
Microstresses 

In the following V(t) C Q shall denote an arbitrary volume that is transported with 
the flow. 

In order to describe the change of the free energy due to diffusion, we introduce 
a chemical potential p and the outer unit normal v to dV{t) such that 




pj-vds. 



is the energy transported into V(t) by diffusion and ds x denotes integration with 
respect to the surface measure. 

Surface forces: Moreover, we assume the existence of a generalized (vectorial) surface 
force £ such that 




<pi-uds.. 



represents the working due to microscopic stresses. Above and in the following tp is 
the material derivative d t p + v • Vp. Finally, we note that 

/ (Tv)-vds x 

JdV(t) 
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describes the working in a given volume V(t) due to the macroscopic stresses in the 
fluid. 

Second law of thermodynamics/local dissipation inequality: As in [20J, we assume the 
following dissipation inequality, which is the appropriate formulation of the second 
law of thermodynamics in an isothermal situation: 

d 
dt 



e(v, (p, V<£>) dx < / (Tv)-vds x + / (p£-vds x — / p3-vds 2 

V(t) ' JdV(t) JdV(t) JdV(t) 



for every volume V(t) transported with the flow. This means that the change of 
total energy in time is bounded by the working due to macroscopic and microscopic 
stresses and the change of energy due to diffusion. 

We recall the transport theorem, see e.g. Liu [2U Theorem 2.1.]: 

at JV(t) JV(t) JdVit) JV(t) 

where the exterior normal velocity of V(t) is given by v ■ v(t) on dV(t), i.e., V(t) is 
transported with the flow described by v. Therefore the equivalent local form is 

d t e + v • Ve - div(T • v) - div(^) + div(/iJ) =: -V < 0. (2.9) 

Using ffTTTU]) . f !23|) and we will simplify V. First of all, multiplying ([ZH]) with 

v and using (12 .7p . one easily obtains 

dt ( P ^2~) + V ' V ( P ^2~) = div(T ' V) ~ T : W " ^~2~ (2 ' 10) 

since 

( |v| 2 \ / |v| 2 \ Ivl 2 

d t (pv) ■ v + div(pv ® v) • v = d t [p— ) + v ■ V I p— j + {d t p + v • Vp) — . 

Moreover, 

a , / + v .v/ = ^ + ^.(v P y (2.11) 

and 

<9pc 

div(0£) — div(/iJ) = div £ + V(<i>) • $, + yU— — — V/x • J 

= (||// + div £V + (W • £ + Vv : (Vip <g> - V/i • J, 
because of (12.71) and where we have used 

d Xj (p = d t d Xj ip + v • Vd Xj ip + d Xj v ■ Vy? = (d x .tp)' + d Xj v ■ Vy?. (2.12) 

9 



Thus we conclude that (12.91) is equivalent to 

\<9v? <9v? dtp 2 J 

+ ( Jf- - ^] (Vv?)' - (S + Vv? <g> : Vv + V// • J < 0, (2.13) 



where we have used divv = 0. 

In order to motivate the constitutive assumptions, we will derive some restric- 
tions for the constitutive relations specifying S,/i, J,£ by an argument typical for 
rational continuum mechanics: To this end, we assume that S, J,£ are functions of 
Dv, ip, Vv?, /i, V/i only. Moreover, we assume that S is symmetric in order to have 
conservation of angular momentum. Invoking general external forces and mass sup- 
plies in the equations, one argues that <p, <p, Vv?, (Vv?)', /i, V/i, (Vv — |divvl) can 
attain arbitrary values for a given point in space and time and since / and £ do not 
depend on ip, (Vv?)', we conclude from (I2.13P that 

necessarily. In particular, £ depends only on <p, Vv?. Hence (12 . 13[) reduces to 

-(S + Vv? ® £) : Dw - (Vv? <8>£) : ^(Vv- Vv T ) + V/i • J < 0, (2.14) 

where Dv = |(Vv + Vv T ). Since the skew part of Vv can attain arbitrary values 
independent of Dv, we conclude that 

£ ® Vv? = Vv? ® £. 

Thus |£| 2 |Vv?| 2 = (Vv? • £) 2 and therefore 

df 

£(<p, Vv?) = a(v?, Vv?)Vv? = g^(<P, Vv?) (2.15) 

for some a(ip, Vv?)- 

On the other hand, the first term after the equality sign in (I2.14p is linear in (p 
and /, J and /i are assumed to be independent of (p. Therefore the first term in (12.141) 
has to vanish to satisfy (I2.14p in general and 

dpc df dp I v 1 2 

V = di ~ divK ^' " d^~ ■ 

Finally, the local dissipation inequality is satisfied if and only if 

V= (s + Vv?®t!^- J : £>v- V/i- J > 0. (2.16) 

V <9Vv?y 
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Here S + V<p ® is also called viscous stress tensor since it corresponds to irre- 
versible changes of the energy due to friction in the fluids. 

Constitutive assumptions: Motivated by Newton's rheological law, we assume that 

df 

S + V<p ® ^=f - = 2r){cp )Dv 

for some function f)((p) > 0. 

Finally, we choose J(<p, V<p, p, Vp) in the form 

J(<p, Vp, p, Vp) = -m(<p)Vp, 

where m(<p) > 0, which corresponds to a generalized Fick's law. We remark that J 
can be chosen to be nonlinear with respect to Vp as long as (I2.16P is fulfilled. 
Summing up, we derived the following diffuse interface model: 

d t {pv) + div(pv <g> v) - div (2r)((p)Dv) + Vp = - div(a(<p, V<p) Vp ® V<p), (2.17) 

divv = 0, (2.18) 
dtpc(<p) + v • Vpc(<p) = div(m(<p)Vp) (2.19) 



together with 



|^p = ^(<p, Vp) - div(a(p, Vp)Vp) - ^ W!, (2.20) 



where a(<p, V<p) satisfies (12.15p . Assuming the normal component of v and the normal 
derivative of p to vanish on the boundary of the fluid domain Q, f n pc(-,t) dx is 
a constant in time. By (12. 4p and (12. 5p . we may express both mass densities as 
affine linear functions of pc. Hence, the total mass J n pidx, i = 1,2, of each liquid 
component is conserved. 

If p = c is the mass concentration difference, then pc(c) = p(c)c, and we obtain 

<9t(pv) + div(pv <g> v) - div (2r)(c)Dv) + Vp = - div(a(c, Vc) Vc ® Vc), (2.21) 

divv = 0, (2.22) 
<%(pc) + v • V(pc) = div(m(c)Vp) (2.23) 



together with 



d(pc)^ = df^ ^ _ diy(a(C) Vc)Vc) _ |Phf (2>24) 



(9c (9c 9c 2 

Here 

<9(pc) 2 1 1 

— — = ap with a = — + — . 
oc 2pi 2p 2 
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In the case that (p — pc is the density difference, we have pc(p) — (p, |p — 1, and 
therefore 

d t (pv) + div(pv ® v) - div (2r](p)Dv) + Vp = - dw{a{<p, Vp)Vp <g> Vtp), (2.25) 

divv = 0, (2.26) 
d t <p + v • Vv? = div(m(</?)V/i) (2.27) 



and 



|2 



df do Ivl 

// = -±(<p, V<p) - div(a(<p, V<p) Vy?) (2-28) 

Finally, in the case that p = m 2 — U\ is the difference of volume fractions, we have 
pc(cp) = + P±±Picp and we obtain the system (l2~25j) . together with 

%> + v- v> = -divJ v , (2.29) 

where we use a rescaled flux 3 V = ( /?1 ^ 2 ) 1 J, and a rescaled chemical potential 
Hip = pl ^ p2 pL. We hence obtain 

J v = -m(tp)Vp, v , 

where m(<p) = ( Pl ^ P2 ) fh(<p). Usually we take the difference of volume fractions as 
order parameter. This has the advantage that the mass difference depends linearly on 
the order parameter and as usual in phase field models the values <p = ±1 correspond 
to unmixed "pure" phases. 

2.2 Derivation based on a Local Dissipation Inequality and 
the Lagrange Multiplier Approach 

It is also possible to exploit a dissipation inequality without introducing generalized 
surface forces. We now assume a dissipation inequality 



4 / e(v,(p,V(p)dx+ [ 

at JV(t) JdV(t) 



J e • v ds x < 



for every volume V(t) transported with the flow for some general energy flux J e , 
which will be specified later. 

Then the equivalent local form is 

<9 t e + v- Ve + divJ e < 0. (2.31) 

Because of the conservation law ( I2.3p . we conclude that for every scalar function 
the inequality 

- V := d t e + v • Ve + div J e - \ {d t {pc) + v • V(pc) + div J) < (2.32) 
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has to be valid. 

Using (EE}, (PI]) , flZZEU) we obtain that f l2~32j) is equivalent to 
/<9/ 9/) |v| 2 <9(pc) \ . 

+ ( (V^)" + div (T • v - A V J - J e ) — S : Vv + VA^ • J < 0, 

where we have used divv = 0. Making use of (12.121) we conclude that the latter 
inequality is equivalent to 

+ div ( T ■ v - A V J + - J e ] - ( S + V<^ g> -^J— ] : Vv + VA^ • J < 0, 



where again the equation div v = was used. 

If we assume now that S, A^,, and T ■ v — A^J + g^y? — Je are independent of (p, 
we conclude that the first term after the equality sign in (I2.33P has to vanish for all 
values of (p. Thus 

,d{p± = d£ _ f df \ 

* dip d<p dtp 2 V \dVip ) ' 

If we denote A^ = /i, then we obtain the same identity for the "chemical potential" 
fi as before. Moreover, if we now specify the energy flux as 

J e = T.v-^J + ^, 

we end up with the local dissipation inequality 

V = ( S + ® J : Dv — V/i • J > 0, 

V dVipJ 

which is just (I2.16p . Thus we can derive the model (I2.17p -( l2~20|) as before by making 
the same constitutive assumptions. 

3 Onsager's Variational Principle — a Third Ap- 
proach to Derive Thermo dynamically Consis- 
tent Models 

In this section, we follow a third pathway to obtain diffuse interface models in the 
spirit of (I2.25p -( l2.27p and in agreement with the postulations of thermodynamics. It 
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is based on Onsager's variational principle, see [27] and [TT], [IB], [28J for applications 
in multi-phase flow. To widen the range of applications, we discuss two additional 
features. We allow for gravitational forces or, alternatively, we include the transport 
of a soluble species across fluidic interfaces as an additional effect. For simplicity, we 
refrain ourselves in the second case to a species that does not influence the surface 
tension at the interface. Adopting the notation from the previous sections, the order 
parameter will be denoted by <p and we take it to be the difference 1*2 — 1*1, of the 
volume fractions Uj, j = 1, 2, of the two liquids involved. Hence, 

PVP) = — 2 — V + — 2 — ' ' 

where p\ and p2 are the specific densities of liquid 1 and 2, respectively. With w, we 
denote the mass density of the soluble species which we assume to be dilute. 

Our starting point is the following set of general evolution equations, see Section [2J 

<9t(pv) + div(pv <g> v) - div S + Vp = K, (3.2) 

div v = 0, (3.3) 

%> + v • V<p + div J v = 0, (3.4) 

d t w + v ■ Ww + div J w = (3.5) 

with p(-) as in (13. ip and J v = ( pl t p2 ) 1 J is a rescaled flux, compare the discus- 
sion after (I2.28p . The stress tensor S is symmetric (due to conservation of angular 
momentum), and K denotes the force density. The additional equation (13.51) is the 
mass balance of the soluble species and J w is the corresponding mass flux. These 
equations are supposed to hold in a space-time cylinder fl x (0, T) with Q C M. d being 
the domain where the process takes place. Conservation of mass requires that the 
normal components of and of 3 W vanish on dQ. As free energy, we choose 

f(<p,V<p)dx+ / {g(w) + (3(ip)w}dx 

with f(z,p) = fi(z) + feip)- Here, g(-) is an entropic term, and P((p) attains for 
(p < — 1 or <p > 1 the values (3i or /?2, respectively. In a sharp- interface limit, these 
parameters will reappear through the Henry jump condition 

^ = exp(/3 2 - A) (3.6) 

at the interface separating the two phases (cf. Subsection I4.3.4p . The total energy is 
given by the sum of kinetic and free energy, hence 

F= j PM| v |2 rfx+ j f^y^ dx+ f {g( W ) + p( V ) W }dx. 

Jn * Jn Jn 
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The time derivative of the free energy is given as 

v\ 2 d t <pdx + J ^-d t ipdx- J div (^^ d Mdx + J p(v)v ■ d t v dx 
+ I (g'{w)+ p(ip))d t wdx+ / (3'(ip)wd t ipdx, 



dF _1 r dp 
dt 2 J n dip 



where we assumed that the normal part of -J^ vanishes on dVt. Observe that 

if d 
d t v = - < K + div S — -rr-dM v — divfpv <g> v) — Vp 
p { dip 

and 

/' /' I 

i dx 



f v • div(pv <S> v) dx = — ( -|v| 2 v-Vp< 
Jn Jn 2 



if v = on <9f2. The last identity follows from integration by parts using divv = 0. 
Inserting the evolution equations and integrating by parts gives, assuming v = on 

on. 



dt 



1 I l^l v | 2 ( v ' Vy? + div 3 v )dx- [ |^(v- Vy? + div J ) dx 

2 Jn * J n d<p 

+ J div (y^^j (v • Vy? + div 3 V ) dx - j f3'(<p)w (v • Vy? + div J v ) dx 

+ / v • { K + div S — TT-d t ip v — divfpv ® v) — Vpl dx 
Jn I d<p ' J 

- J (g'(w) + /%)) (v-Vw + div J w ) rfx 

+ f v-Kdx-]- [ v - V(p^f-\v\ 2 dx - [ v-Vw (^(w) + /%)) da; 
in 2 J n dip J Q 

+ / 3 V ■ Vfx v dx — Dv : S dx + / J«, • Vp„, dx. 
Jn Jn Jn 

Here, we used 

- / |v| 2 |^%>cfe = / (v Vy? + div Jj|^|v| 2 cfe 
and we abbreviated 

f df \ Of I dp 

and 

p w := g'(w)+p(<p) 
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to denote the chemical potentials corresponding to (p and w, respectively. Recall 
that in general entropy production is due to external force fields and to gradients 
of velocity and of chemical potentials, see [8], Chapter 3. However, if the specific 
densities of the external forces^] acting on the different species are identical, then 
those forces do not contribute to entropy production. Therefore, if no other external 
forces than gravity forces are applied, we may identify the rate of change of the 
mechanical work with 

dW f _ , f _ , f 1 f _ dp 



dt 
As 



/ v-Vpp^dx — ( v-Vwfi w dx+ ( v-~Kdx — — j v- Vp-^-\v\ 2 dx. (3.7) 
Jn Jn Jn 2 J Q dip 



J In K-grav ■ v dx case with gravitational forces 
dt 1 case with soluble species, 

equation f)3.7p is satisfied if 

{p^Vp + !f^l v | 2 Vy3 + K grav case with gravitational forces, 
H^Vp + fi w Vw + !f^l v | 2 Vy3 case with soluble species. 

Here, K grav denotes the gravitational force. 

To determine the fluxes J^,, 3 W and the stress tensor S, we introduce the dissipation 
functional 

$(J V ,J^,S):= f IJ^L + ^L + J^-ldx. 



\ 2M((p) Ar](ip) 2K((p) 



w 



We use Onsager's variational principle which postulates 

<*(j„,j„,,s) (^~( J ^> J -> S ) + $ ( J ^' J -' S )) = °' 

This gives 3 V = -M((p)Vfi v , J w = -K(p)wVp w , and S = 2r)((p) Vv + 2 VvT . Alto- 
gether, when a soluble species is around, we end up with 

d t (pv) + div(pv <g> v) 

-div(2n((p)Dv) +Vp = ^Vp + ^ w Vw + ]~\w\ 2 Vp, (3.8) 

divv = 0, (3.9) 
d t <p + v- V(p-div(M(<p)Vp: v ) = 0, (3.10) 
d t w + v • Vw - div (K(ip)wVfi w ) = 0, (3.11) 

p w = g'(w)+/3(ip) (3.13) 



2 i.e. the quotient of force density and mass density 
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where Dv = |(Vv + Vv T ) and p{tp) = P2 2 pi (p + Pl + p2 . The evolution equation for 
w becomes 

d t w + v - Vw = div{K((p)wV{g'{w) + /%))). (3.14) 
In the case of gravitational forces, the system has to be changed accordingly. 
Remarks: 

• If we choose g(w) = w(\ogw — 1), the equations (|3.1ip .( )3~7l3"j) result in the 
diffusion equation 

d t w + v • Vw = div(K(ip)(Vw + ioV/%))). 

• The interfacial force term K = /i^Vy? + fi w Vw + |f^l v | 2 Vy5 can equivalently 
be written as 

V(f(<p, V<p) + g(w) + p(tp)w) - div (\?<p (g- Dj{ ~'-- 



Here, the first term is of pressure type whereas Vy?®^^ provides an additional 
stress tensor contribution representing interfacial forces. We can hence conclude 
that the derivations in Sections [2] and [3] up to a reinterpretation of the pressure 
lead to the same diffuse interface model. The analogous observation for "Model 
H" has been already discussed in [20J. 

It is also possible to derive (I3.8p -( 13~T3|) with the approaches discussed in Section 
® 

Taking the definitions of // v and K into account, we observe that K in fact 
does not depend on v. 



4 Sharp interface asymptotics 

In this section we identify the sharp interface limit of the diffuse interface model 
introduced in the preceding sections. We will use the method of formally matched 
asymptotic expansions where asymptotic expansions in bulk regions have to match 
with expansions in interfacial regions. There are four different asymptotic limits of 
interest. Two use a constant mobility and will either lead to a model where diffusion 
takes place through the bulk or to a model without any diffusion through the bulk, 
see also Abels and Roger [3] for the case when the densities in the two fluids are the 
same. Two cases are based on a mobility which is zero when the phase field takes 
the values ±1. In this case one of course does not observe diffusion through the bulk 
in the sharp interface limit. But depending on the scaling we will either see surface 
diffusion along the interface or not. 
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4.1 The governing equations 

As usual for phase field models we introduce a scaling for / with respect to a small 
length scale parameter e as follows 

/(y?,Vy.) = ^|Vyf + -^) 
2 e 

where a is a constant related to the surface energy density. As in Section 3, we choose 
the difference of volume fractions as order parameter and we consider the following 
system 

d t (pv) + div(pv ® v) - div(2r? (</?)£> v) + V> = -be div( Vy? ® V<p) , (4.1) 

divv = 0, (4.2) 

d t ip + v • V<p = div(m £ (<p) V/i) , (4.3) 

/i = " oe^p - p'Mllvl 2 + ^(^)» , (4.4) 

d t w + v • Vw = div{K(ip)wV(g'(w) + /%))). (4.5) 

To simplify the notation we drop the ip as index in the chemical potential. We assume 
that 

. ^) = ^+^, 

• g is convex, 

• p(tp) is smooth with 0(1) = /3 2 , = 

• ?y(<p) is smooth and positive with 77(1) = r] 2 , = rji, 

• ip((p) is a double- well potential such that = 1) = and ^(z) > if 
z <£{!,-!}, 

• K(ip) is smooth such that K(l) = K 2 , K(-l) = K x . 
For the mobility m £ we distinguish four cases: 



m £ (ip) 



where m ,mi > are constants and (.)+ is the positive part of the quantity in the 
brackets. The total relevant energy in this scaling is 

F e (<p, v, w) = J |^|v| 2 + y I Vy?| 2 + -^W) + g(w) + 0(<p)w} dx . 
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For a solution (v e , p £ , <p £ , fi e , w £ ) of the system fl4.1|) - fl4.5p we perform formally matched 
asymptotic expansions. It will turn out that the phase field ip £ will change its values 
rapidly on a length scale proportional to e. For additional information on asymptotic 
expansions for phase field equations we refer to [151 E] • 

4.2 Outer expansions 

We first expand the solution in outer regions away from the interface. We assume an 

expansion of the form v e = Y^T=o £kv k, ip e = YlT=o £k{ P k ^ ^n expansion of f)4.4p 

in outer regions gives to leading order i[>'(ipo) = and we obtain the stable solutions 
±1. We will denote by Q ± the regions where ipo = ±1. The leading order expansion 
of the other equations are straightforward. We obtain: 



d t w + v ■ Vw = KiV ■ (w Vg'(w )) , (4.9) 

where i = 1,2 for x G and p\ = p(— 1) = pi,p2 = = p 2 . Due to the 

divergence free velocity we obtain 2 div Dv = Av and hence (14.61) simplifies to 

Pi(«9<v + div(v <g> v )) - ?7iAv + Vp = . 

We remark that (14. 9p leads to the convection diffusion equation 



in the case that g(w) = u;(log w — 1). In the cases II-IV we will not need the chemical 
potential p in the bulk. 

4.3 Inner expansions 

We now make an expansion in an interfacial region where a transition between two 
phases takes place. 

4.3.1 New coordinates in the inner region 

Denoting by T = (T(t)) t > the smoothly evolving interface which we expect to be 
obtained in the limit when e tends to zero, we now introduce new coordinates in a 
neighborhood of T. Choosing a time interval IcK and a spatial parameter domain 
U C R d_1 we define a local parameterization 



of T. By i> we denote the unit normal to T(t) pointing into phase 2 (which is the 
phase related to <p = 1). Close to j(I x U) we consider the signed distance function 



Pi(d t Vo + div(v <S> v )) - 2r)i div Dv + Vp 

div v 

in case I Ap 



0. 

0. 



(4.6) 
(4.7) 
(4.8) 



d t w + v • Vw = KiAw 



j:IxU^R d 
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d(t,x) of a point x to T°(t) with d(t,x) > if x G f2+(t). We now introduce a local 
parameterization of / x R d close to 7(1 x [/) using the rescaled distance z = - as 
follows 

G £ (t, s, z) := (t, j(t, s) + ezu(t, s)) . 

We denote by 

V = Off - v 

the (scalar) normal velocity and observe that the inverse function (t, s, z)(t,x) := 
(G e ) _1 (t, x) fulfills 

To derive the last identity we used (2.6) and (2.20) of [7\. For a scalar function b(t, x) 
we obtain for b defined in the new coordinates via b(t, s(t,x), z(t,x)) = b(t,x) the 
identity 

^-b(t, x) = d t zd z b + d t s ■ W s b + d t b = ~Vd z b + h.o.t. 

where h.o.t. stands for higher order terms. With respect to the spatial variables we 
obtain, see Appendix, 

w x b = w T J+\d z bu (4.10) 

where Vr ez is the tangential gradient on 

T £Z := { 7 (s) + ezu(s) | s G U} 

where here and in what follows we often omit the t-dependence. For a vector quantity 
j(t, x) written in the new coordinates via j(t, s(t, x), z(t, x)) = j(t, x) we obtain 

V ie -j = Vr„-j + i9J-i/, (4.11) 

where Vr ez • j is the divergence of j on T £Z . In the Appendix we compute 

A x b = A r J - i(« + sz\S\ 2 )d z b + ±d zz b + h.o.t. , (4.12) 

where k is the mean curvature (the sum of the principal curvatures) and |«S| is the 
spectral norm of the Weingarten map S. In addition we note that (see Appendix) 

Vr E j(s,2) = V r &(s,z) + h.o.t. , 
Vr ez -j(s,2;) = V r -j(s,z)+ h.o.t. , 
A r£z b(s,z) = A r S(s, z) + h.o.t. 

where Vr, Vr-, Ar are the surface gradient, the surface divergence and the surface 
Laplacian on T. 



20 



4.3.2 Matching conditions 

We now assume an e-series approximation of the unknown functions ip, fi, v,p, w, . . . 
which in the inner variables we will denote by M, V, P, W, . . . . Denoting by $o + 
e$i + . . . the inner expansion and by ip + eipi + . . . the outer expansion of the phase 
field we obtain the following matching conditions at x — "y(s): 

lim $oO,s) = <p (x±) , (4.13) 

z— y±co 

lim d z $x(z,s) = V<po{x±)-v (4.14) 

z— >±oo 

for z — > ±oo, where (po(x±), • • • denotes the limit lim<£>o(^ ± Su). In addition we 

<5\0 

obtain that if $1(2, s) = A±(s) + B±(s)z + o(l) as z — > ±00 the identities 

A±(s) = <pi(x±), B ± (s) = V<p (x±) ■ v (4.15) 

have to hold (see |14j . [18]). Of course similar relations hold for the other functions 
like v, //, 

4.3.3 Deriving the sharp interface limit 

Plugging the asymptotic expansions into fl4.ip -( )4~5]) we ask that each individual co- 
efficient of a power in e vanishes. The equation (14.41) gives to leading order ~: 

O = c^$ -V'($o) (4.16) 

and from (14. 13j) we obtain 

$o(z)->-±l for z^ioo. (4.17) 
We now choose the unique solution of (I4.16p . (I4.17P which fulfills 

$o(0) = 0. 

We in particular obtain that $0 does not depend on t and s. Equation (14. 2p gives to 
leading order 

d z V -v = d z (V -v) = 0. (4.18) 
The matching condition requires that (Vo • v){z) is bounded. Hence 

(v -i/ )(x+) = lim (V • u )(z) = lim (V • u )(z) = (v • u )(x-) . 

z— >oo z— 00 

This implies 

[v o -u]t = 0, 

where [u]l(x) = u(x+) — u(x-) denotes the jump of a quantity at the interface. 
Applying ( I4.10p for each component we obtain 

V x v = i^V<8»i/ + Vr„V, 

D x v = H(9*V®i/ + i/®9 z V) + i(Vr e ,V + (Vr„V) T ). 
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Defining £ (A) = |(A + A T ) for a quadratic matrix A we compute 

V x -{n(<p)D x v) = ^(Ty^f^V® + i^(77($)f(Vr 8 ,V))i/ 

+ ^V r „ • (v(®)£(d z V ® i/)) + Vr £Z • (iK*)£(V r «,V)) 
= 4^(?7($)£(^V ® + i&(77($)£(V r „V)i/) 

+^V r „ ■ fa(*)£(a,V ® i/)) + Vr £2 • (t/(S)5 (V r „ V)) , 

where we used <9 2 z/ = 0. We conclude from (14. 18j) 

(i/®^V )i/=(^Vo-i/)i/ = 0. 

The fact that $o does not depend on t and s and f)4.10p imply 

V<p <g> V<p = ^(d z $ )V (8) i/) + i<9 z $id z $ (> ® v) + h.o.t. 

and since (Vrf) ■ ^ = we get 

£V ■ (Vy? (8) V<p) = ^(d&ofv + i(d 2 $ ) 2 (Vr ■ + \d z {d z ^d z ^ + h.o.t.. 
Hence we obtain at the order from the momentum equation 

&d t (d t ^o) 2 v + dM$ )d z V ) + (d z P„i)v = . (4.19) 
Multiplying (I4.19P with u, taking d z u = and d z V ■ v = into account gives 

&d z ((d z <5> ) 2 ) + d z P^ = . 

Hence (I4.19P implies 

dM$o)d z V ) = 0. (4.20) 

Matching implies that V (z) is bounded. This implies that (14.201) interpreted as an 
ODE in z has only solutions Vo which are constant in z. This implies after matching 

[v ]t = 0. 

We now want to analyze the momentum equation to the next order. The term 
V • (r)(<p)Dv) gives to the order -, 

d z (r ] (%)S(d z V 1 ® u)u) + dM<S> )£(V r V )v) . 

Matching requires lim d z \~i(z) = Vv (x±)y and hence 

z— >±oo 

d z \i ®v + V r V -»■ V x v for z -»■ ±oo . (4.21) 
Altogether we obtain for the momentum equation at order -: 

- ^(p($o)Vo)V + 9 2 (p($o)(Vo®Vo))^ 

- 2d z (r ] (%)£(d z V 1 ® - 29 z (r / ($ )^(V r V ) l /) 

- a{d z %) 2 Ku + <7^(d z $ 1 <9 z $o)f + <9 Z P ^ = . 
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Integrating with respect to z gives after matching and using (I4.2ip 

-[Pov ]lV + [p v ]^v • v - 2fo£(V x Vo)]ii/ - a (J (d z %) 2 dz\ kv + \po]±v = . 

Defining c := (d z $> ) 2 dz and a = c <j we obtain 

[pov ]± (v • v - V) - 2[riDv ]tv + [po]±> = okv . 
4.3.4 The diffusion equations to leading order 

The analysis of (14. 3 p now depends on the ansatz for the mobility. We have to distin- 
guish between four cases. 

Case I : m E ((p) = m . 

At order -4 we obtain from (14.31) 



= d z (m d z M v) ■ v = m d zz M . 

Matching implies that M is bounded and hence M is constant. In particular we 
derive 

Case II : m e (<p) = em . 

At order | we conclude from ( 14. 3 p 

- Vd z % + (v • v)d z <5> = d z {m d z M u) ■ v = m d zz M . (4.22) 

Since = 0, we obtain from matching 

d z M ->■ for z -»■ ±oo . 

Integrating (I4.22p with respect to z gives 

V = \q ■ v . 

In addition we get d zz M = and hence M does not depend on z. 

Case III : m e ((p) — ra\{l — y? 2 )+ . 
At order ^ we get 

= d z ( mi (l - <5>l)d z M v) ■ v = d z {m x {\ - &*)d z M ) . 

Matching implies 

mi(l - ®l)d z M ^ for z ->■ ±oo . 

We hence obtain 

7^(1-^)9^0 = 0, 
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which implies 

M = M (s,t). 

Case IV : m e ({p) = emi(l — y? 2 )+ ■ 
At order | we get 

-V<9 2 $ + (V ) • ud t $ = d z ( mi (l - & )d z M ) 

and hence combining arguments from the Cases II and III above we obtain 

V = v • v and M = M (s, *) . 

We now analyze the diffusion equation for the soluble species. Equation (14. 5 p gives 
to leading order \ 

d z {K{$ )Wod z {g'{Wo + /3($o))) = . (4.23) 

Matching to the outer solution gives that g'(Wo) + /3($o) is bounded. Multiplying 
(I4.23P by g'(Wo) + /3($o), integration and integration by parts gives 

K(%)W \d z (g'(W ) + (3(<5> ))\ 2 dz = . 

Assuming Wq > we obtain that 

g'(W )+P($ ) 

does not depend on z. Hence 

[^H) + /%o)]i = o. 

This implies 

logw (a;+) - logiwo(ac-) + P2 ~ Pi = 
which yields the Henry jump condition 

wpjx-) 

— — — = exp(/3 2 - Pi) • 
w (x+) 

4.3.5 The generalized Gibbs— Thomson equation 

The equation for the chemical potential gives to the order e° 

M = aV"($o)$i - od zz ^ x + &d,$ K - ip'($ )|V | 2 + P'(%)W . (4.24) 

In order to be able to obtain a solution $1 from (I4.24p a solvability condition has to 
hold. This solvability condition will yield the generalized Gibbs-Thomson equation. 
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We multiply (I4.24p with <9 2 $ , integrate with respect to z and obtain (using the facts 
that M and Vo do not depend on z): 

/oo 
(^'($o)^$ <f i - d zz ^d z %)dz 
-OO 

/oo p oo /*oo 

(<9 2 $ ) 2 ^ - ||V | 2 / d zP {%)dz + / d z /3{%)W dz. 
-oo J —oo J —oo 

We obtain after integration by parts, using the fact that d z §o(z), d zz §o(z) decay 
exponentially for \z\ — > oo, 

/oo 
d g (i//($ ) - d zz $ )$idz + ac K - i|V | 2 [po]± 
-OO 

poo POO 

+ / d x (p($ )W + g(W ))dz- / (P($o)+g , (W ))d g W dz. 



Since /3($o) + {/(Wo) does not depend on z and since ip'($o) — <9 22 <3>o = 0, see ( I4.16p . 
we obtain 

2/i = ok - ||v | 2 [p ]t + [g(w ) - g'(wo)w ]t ■ 

This identity has been derived for a general function g. In the case g(w) = w log w — w 
we obtain 

2/i = crre - ||v | 2 [po]- - [wo]l • 

4.3.6 Interfacial flux balance in the sharp interface limit 

We now expand the equations (14. 3 p and (14. 5 p further in order to obtain contributions 
of the diffusive fluxes in the interface. The result will depend on the choice of the 
mobility. In the cases II and IV the interfacial diffusive fluxes for <p are scaled such 
that they do not contribute to a limiting sharp interface problem. We hence consider 
only the cases I and III. 

Case I : m s ((p) = m . 

At order | we deduce from (14. 3p 

(-V + V ■ u)d z % = d^modM) . (4.25) 

Matching gives d z M 1 — > Vpo ■ v for z —> ±oo. Integrating ( I4.25P gives 

2(-V + v -i>) =m [V// -f]i- 
Case III : m e (<p) = ^(1 - <p 2 ) + . 

We have up to order e° in the interfacial region (setting m((p) = TOi(l — if 2 )) 

V ■ (m(ip)Vfi) = ±d z (m($ )d z M ) + 

id^m^dM + m'($ )$AM ) + ±V r • (m($ )d z M v) 

+d z {m{%)d z M 2 + m'^Q&Mx) + V r • {m{$ Q )d z M x u) + V r • (m($ )V r M ) . 
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Using d z M = we obtain from (14. 3 p at order 4^ 

= 9 z (m($ )9,M!) . 

Matching gives 

m($o)d z Mi -> for z — > ±oo . 
Hence m( <£>o)<9 2 Mi = and 

M x = Mi(s,t) . 

At order ~ we obtain from (I4.3p . using d z M\ = and d z Mo = 0, 

- Vd z % + (V ■ u)8M = d z (m(%)d z M 2 ) + V r • (m($o)V r M ) . (4.26) 
Matching gives 

m($ )d z M 2 -> for z -> ±00 . 

Integrating (I4.26P gives 

/oo 
Vr-(m($o)VrM )dz. 
-00 

Since m(<& ) does not depend on s, we obtain Vr • (m($ )VrM ) = m($ )^r^o- 
Altogether we deduced 

2(— V + v ■ 1/) = mArfjio 

where rh = m($>o)dz. 

Finally we deduce the flux balance for w at the interface in the limit e — > 0. 
Equation (14.51) gives to order | 

- Vd z W + (V • f)<9 z Wo = d s (K($o)W d z (g"(W )Wi + /3'($o)$i)) • (4.27) 
Matching gives 

d z (g"(W )W 1 + /3'($o)$i) ->• VQ?>o) + /3(±1)) ■ i> 
for z — >■ ±00. Integrating (I4.27P gives 

(V - v • 1/) Hi = -{Kw Vg'(w )]t ■ u. 
In the case g(w) = wlogw — w this identity reduces to 

(V-v -u)[w Q ]t = -[KVw ]t-u. 
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5 Free energy inequalities for the sharp interface 
limit 



The sharp interface limit is now given as follows. We search for a smoothly evolving 
hypersurface T which for all t > separates Q into open sets fi_(t) and and 
we look for functions v,p,w which are all defined for x E Q and t > 0. In the cases 
I and III we seek in addition a chemical potential \i which in case I is defined on 
Q x (0, oo) and for case III the potential fi is defined on the interface T only. 

Assuming for g the form g{w) = w(\ogw — 1) we need to solve in all four cases 
the system 

Pi(d t v + div(v <g> v)) - t^Av + Vp = 0, 

divv = 0, 
d t w + v • Vu> = KiAw 

in the bulk regions f2_ and Q + . In case I we in addition solve 

Afj, = 

in the bulk. On the interface T we require 

0, 

exp(/3i -&), 
-[KVw\t ■ v 

where w 1 = w(x—) and w 2 = w(x+). All other conditions depend on which of the 
cases I-IV we consider. 

Case I (m £ (ip) = m ): 

[p]lv(v • is - V) - [2rjDv]±v + \p]±v = <tkv, (5.1) 

2(vi/-V) = 771o[Vax]1 • i/, (5.2) 
2^ = o-k - ||v| 2 [p]l - (5.3) 

Remarkable here is that the additional term [p]lv(v- v — V) enters the stress balance 
at the interface. We also note that the non-standard term ||v| 2 [p|l enters the Gibbs- 
Thomson law. 

Cases II and IV (lim^o m s ((f>) = 0): 

— [2rjDv]t • v + = gkv , 

V = vi/, 

i.e., in this case we recover a standard free boundary problem for the Navier-Stokes 
system which in addition is coupled to the flow of a soluble species. 



Ml = 

w 2 

Wi 

(v-v-i/)Ml = 
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Case III (m e ((p) = m\{\ — <f 2 )+): 
Beside (15. II) and (15. 3p the identity 



2(v • v - V) = mA r n 

has to hold. In this case the diffusion of the two components is limited to the 
interfacial region. In fact the well-known surface diffusion flow V = Apft, see [2"Bl El 
[T2"| [T3] . is in the new model coupled to fluid flow. 

In all cases we require v = 0, Vw • v — on dfl and in case I we in addition 
require V/z • v = on dQ. In order to derive a free energy inequality for the sharp 
interface limit we need the following transport identities, for a proof see [7]. 

Lemma 5.1 (Transport identities). Let T = (T(t))t>o with T(t) C £1, for all t > 0, 
be a smooth evolving hypersurface and let f be a quantity which is smooth for all 
t > and x 6 O \ T(t) and such that [/]! exists at T. Then it holds 



dt 

and 



— ( lds x — — f kV ds x — — f ■ (W) ds a 

Jr(t) Jr(t) Jrtt) 



T(t) Jr(t) Jr(t) 



| ; / fdx= [ dtfdx- [ [f]tVds 

at Jn Jn Jr(t) 



We are now in a position to compute the dissipation rate and in conclusion derive 
a free energy inequality. 

Theorem 5.2 (Free energy inequality) . 

A sufficiently smooth solution of the sharp interface problem with T(t) C Q, for 
all t > 0, fulfills 



d 
dt 



/ ( ?l v ! 2 + (g( w ) + Pwj) dx + / ads 
Jo Vz! J Jr(t) 



-V < 



provided the integrals are finite. Here p = p\,p% and (3 = f}\, (3 2 in the two phases. 
The quantity V is given as 

Case I : V = [ 2n| 



Case III : V = / 2rj\ 

Jn Jn w Jr(t) 



\ 2 dx + 


/ ^ 1 


Vw\ 




Jn 




\ 2 dx + 




Vw\ 




In 




\ 2 dx + 


I w 


Vw\ 
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Proof: We now give a proof for the cases I and III in detail and discuss the cases II 
and IV afterwards. For the kinetic energy we compute, using a transport identity, 

d_ 
~dt 



I p^—dx = ! pv-d t vdx- ! [p\t^—Vds x 
Jn * Jn Jv{t) * 

= I v • (— pdiv(v <g> v) + 2?7divDv — Vp) dx — f [p]_-^-V ds x . 
Jn Jv(t) 2 

Using div v = we obtain div(v (g> v) • v = |V|v| 2 • v. Integration by parts on Q + 
and Q- now gives 



37 / P^-dx = - / 27]\Dv\ 2 dx+ [ v • (\\p]t |v| 2 f - [2rjDv)tu + \p]±v) ds x 



dt 



[ \[p\ + M 2 Vds x 
Jv{t) 

i 2rj\Dv\ 2 dx+ / (-|[p]^|v| 2 (v- v- V) + aKV-v)ds 3 
Jn J Tit) 



X*) 

where we used the stress balance (15 .11) and the fact that v is continuous. 

The total free energy of the soluble species fulfills (using g(w) = wlogw — it?) 

d 



, . (g(w) + j3w)dx — / d t w(\og w + (3) dx — / [g(w) + f3w]tV ds x . 
dt Jn Jn Jvit) 

Since d t w + v ■ Vw = KV ■ (wV logw), we obtain 

/ d t w{\og w + /3) dx = / (log w + /?)(— v • Vw + KV ■ (wV \ogw)) dx 
Jn Jn 

= — v • \7(g(w) + f3w) dx — / Kw\ V log w\ 2 dx 
Jn Jn 

— / [(w log w + [5w)KV log w]i • i/ ds x 
Jv(t) 

— — Kw\ V logu;| 2 dx + / [g(u*>) + /3w]lv • f 
Jn Jv(t) 

[{\ogw + /3)-ftTVw]l • Ws* 



— / i^w|Vlogtf| 2 dx + / [<7(io) + (3w]tv ■ uds 3 
Jn Jr(t) 

+ / (logw + /?)(V- v-i/)Hids a 



X*) 

where we used divv = and the fact that Henry's law implies the continuity of 
logw + 0. In addition we have 

— f a ds x = — [ anV ds x . 
dt Jv(t) Jv(t) 
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Altogether we obtain using the generalized Gibbs-Thomson law ( 15. 3 p and the fact 
that g(w) + f3w — (logw + (3)w = —w 



d_ 

dt 



(5.4) 



J (p^TT + 9{w) + /3wj dx + j cr ds x ^j 

= - 2r]\Dv\ 2 dx- / K±\Vw\ 2 dx 
Jn Jn w 

+ [ (vi/- V)(-|[p]l |v| 2 + ok - [w]±) ds x 

Jv(t) 

= - 2r)\Dv\ 2 dx- / K±\Vw\ 2 dx+ / (v • v - V)2// ds x . 
Jn Jn w Jr(t) 

In case I we obtain by integration by parts using (15.21) and A/z = in the bulk 

/ 2(v • v — V)fids x = / m /i[V/i]t • v ds x (5.5) 
Jv(t) Jr(t) 

— — I V ■ (/im V(U) dx = — / m \Vfi\ 2 dx. 
Jn Jn 

In case III we obtain using integration by parts on manifolds without boundary 

/ 2(v • v — V)/x ds x = I m(A r /i)yU ds x = — I m|Vr/-t| 2 ds x . (5.6) 
Jr(t) ' Jv(t) ' Jv(t) 

Combining (l5.4p -f l5T6|) gives the result in the cases I and III. In the cases II and IV 
we have V = v ■ v and the calculations simplify. In particular we do not need to treat 
the term /^^(v ■ v — V)2fids x and diffusion plays no role for the fluid flow. ■ 



Appendix 

We use the notation of Section l4~3l and prove the identities (I4.10p . I4.1ip and (I4.12p . 
Let (si, . . . , Sd-i) G U and G £ x (s, z) = *y(s) + ezu(s), where here and in what follows 
we omit the t-dependence. Then 

<9 S1 7 + ezd Sl v, <9 Sd _ x 7 + ezd Si _ x v, eu 

is a basis of M d locally around V . We define the metric tensor in the new coordinates 
as follows 

9%j = (d Si *y + ezd Si v) ■ (d Sj 7 + ezd Sj v) i, j = 1, . . . , d - 1 , 
9id = 9di = {d Si l + ezds-v) ■ eu = i = 1, . . . , d - 1 , 
g dd = ev ■ eu = e 2 
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where we used that d Si v ■ u = ^d Si \is\ 2 = 0. We set 

Qez = (gij)i,j=l,...,d, Qez = (gij)i,j=l,...,d-l, (Qez)' 1 = (g lj )i,j=i,...,d,(Gez)~ 1 = (g tj )i,j=i,...,d-i 
and hence 



/ 



o \ 



V 





0---0 e 2 



( 



, (Qe 



J 



\ 



V 



o 

0---0 e- 2 



Denoting by Sd the ^-variable we have for a scalar function b(t, x) = b(t, s(t, x), z(t, x)) 



d d-l 

Vsfe = Y,9 l] ds l bd S] G e x =Y,9 l3 d s fid S] G e x + j I d z bd z G\ 

i,j=l i,j=l 

= V v J + l d z bv 



where we used the fact that = for i = 1, . . . ,d — 1. Here Vr £Z & is the tangential 
gradient Vr ez &|r £Z on Y ez := {7(5) + ezv(s) \ s G U}. In addition we compute for a 
vector quantity x) = j(t, s(t, x), z(t, x)) 

d d-l 

V,-j = J29 lJd sJ-d S] Gl=J29 t3 d s J-d Sj Gl + ^d z ]-eu (5.7) 

i,j=l i,j=l 



where Vr ez • j is the divergence on T ez . We remark that Vr ez & • v = , as v is normal 

d z (y T J-v) = 



to r £2 . We hence obtain 



and 



Since <9 z f = 0, we get 



We now compute 



3*(V r „S)-i> + Vr„&.d z i/ = 0. 



Kb = V x ■ (S7 x b) = V x ■ (V r J + \d z b v) 

= Vr £Z • (V r J) + i(Vr«A&) • u 

+ ~d z W Vez ■ v + \(d z S7 r J) • v + ±d zz b ■ v + ±d z bd z v ■ u . 

Because of (Vr £Z d z b) ■ v = 0, (d z Vr £Z b) ■ 1/ — 0, u — V x d, ( 15. 7p and c? z i/ = 0, we 
obtain 

A x 6 = A r J + UA x d)d z b + ±dj . 
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Introducing g £z = det Q £Z we get 

1 d-l 
V9ez . i=1 

Since 

<7y = <9 Si 7 • 9 Sj 7 + h.o.t. , z, j = 1, . . . , d - 1, 

we derive 

V r „S(>,z) = V r 6(s, z) + h.o.t. , 
Vr ez -jO,^) = V r + h.o.t., 

Ar £Z 6(s,2;) = A r 6(s, z) + h.o.t. 

where Vr, Vr- and Ar are computed on T £Z with the metric tensor Q . From Gilbarg 
and Trudinger [19] (Lemma 14.17) we get (denoting by Ki the principal curvatures) 

d-l d-l d-l d-l 

A x c? = > — - = > — = — > ^ — > enjz + h.o.t. 

i=l i=l i=l i=l 

= — k — ez|iS| 2 + h.o.t., 

where k is the mean curvature and |«S| is the spectral norm of the Weingarten map S. 
Hence we obtain 

A x b = A r b - i(« + ez\S\ 2 )d z b + ^9 2Z 6 + h.o.t. . 

References 

[1] H. Abels. Existence of weak solutions for a diffuse interface model for viscous, 
incompressible fluids with general densities. Comm. Math. Phys., 289(l):45-73, 
2009. 

[2] H. Abels. Strong well-posedness of a diffuse interface model for a viscous, quasi- 
incompressible two-phase flow. In preparation. 

[3] H. Abels and M. Roger. Existence of weak solutions for a non-classical sharp 
interface model for a two-phase flow of viscous, incompressible fluids. Ann. Inst. 
H. Pomcare Anal. Non Lineaire, 26(6):2403-2424, 2009. 

[4] F. Boyer. Nonhomogeneous Cahn-Hilliard fluids. Ann. Inst. H. Poincare Anal. 
Non Lineaire, 18(2):225-259, 2001. 

[5] F. Boyer. A theoretical and numerical model for the study of incompressible 
mixture flows. Computers and Fluids, 31:41-68, 2002. 



32 



[6] J.W. Cahn, CM. Elliott, and A. Novick-Cohen. The Cahn-Hilliard equation with 
a concentration dependent mobility: motion by minus the Laplacian of the mean 
curvature. European J. Appl. Math., 7(3):287-301, 1996. 

[7] K. Deckelnick, G. Dziuk, and CM. Elliott. Computation of geometric partial 
differential equations and mean curvature flow. Acta Numerica, 14:139-232, 2005. 

[8] S. R. de Groot and P. Mazur. Non-equilibrium Thermodynamics. North-Holland, 
1969. 

[9] H. Ding, P.D.M. Spelt, and C Shu. Diffuse interface model for incompressible 
two-phase flows with large density ratios. J. Comput. Phys., 226(l):2078-2095, 
2007. 

[10] J. E. Dunn and J. Serrin. On the thermomechanics of interstitial working. Arch. 
Rational Mech. Anal, 88(2):95-133, 1985. 

[11] C. Eck, M. Fontelos, G. Grim, F. Klingbeil, and O. Vantzos. A phase-field model 
for electrowetting. Interfaces and Free Boundaries, 11: 259-290, 2009. 

[12] CM. Elliott and H. Garcke. On the Cahn-Hilliard equation with degenerate 
mobility. SI AM J. Math. Anal, 27(2):404-423, 1996. 

[13] CM. Elliott and H. Garcke. Existence results for diffusive surface motion laws. 
Adv. Math. Sci. Appl, 7(l):467-490, 1997. 

[14] P.C Fife. Dynamics of Internal Layers and Diffusive Interfaces. CBMS-NSF 
Regional Conference Series in Applied Mathematics 53, Soc. Ind. Appl. Math., 
Philadelphia, 1988. 

[15] P.C. Fife and O. Penrose. Interfacial dynamics for thermodynamically consistent 
phase-field models with nonconserved order parameter. EJDE 1995(16): 1-49, 
1995. 

[16] M. Fontelos, G. Grim, and S. Jorres. On a phase-field model for electrowetting 
and other electrokinetic phenomena, submitted for publication. 

[17] H. Garcke, B. Nestler, and B. Stinner. A diffuse interface model for alloys with 
multiple components and phases. SIAM J. Appl. Math. 64(3):775-799, 2004. 

[18] H. Garcke and B. Stinner. Second order phase field asymptotics for multi- 
component systems. Interfaces Free Boundaries 8:131-157, 2006. 

[19] D. Gilbarg and N.S. Trudinger. Elliptic partial differential equations of second 
order. Springer- Verlag, Berlin 2001. 



33 



[20] M. E. Gurtin, D. Polignone, and J. Vinals. Two-phase binary fluids and immis- 
cible fluids described by an order parameter. Math. Models Methods Appl. Sci., 
6(6):815-831, 1996. 

[21] P.C. Hohenberg and B.I. Halperin. Theory of dynamic critical phenomena. Rev. 
Mod. Phys., 49:435-479, 1977. 

[22] H.-G. Lee, J. S. Lowengrub, and J. Goodman. Modeling pinchofT and recon- 
nection in a Hele-Shaw cell. I. The models and their calibration. Phys. Fluids, 
14(2):492-513, 2002. 

[23] H.-G. Lee, J. S. Lowengrub, and J. Goodman. Modeling pinchofT and recon- 
nection in a Hele-Shaw cell. II. Analysis and simulation in the nonlinear regime. 
Phys. Fluids, 14(2):514-545, 2002. 

[24] I-Shih Liu. Continuum mechanics. Advanced Texts in Physics. Springer- Verlag, 
Berlin, 2002. 

[25] J. Lowengrub and L. Truskinovsky. Quasi-incompressible Cahn-Hilliard fluids 
and topological transitions. R. Soc. Lond. Proc. Ser. A Math. Phys. Eng. Sci., 
454:2617-2654, 1998. 

[26] W.W. Mullins. Theory of thermal grooving. J. Appl. Phys., 28:333-339, 1957. 

[27] L. Onsager. Reciprocal relations in irreversible processes I. Phys. Rev., 37:405- 
426, 1932. 

[28] T. Qian, X.P. Wang, P. Sheng. A variational approach to moving contact line 
hydrodynamics. J.Fluid Mech., 564:333-360, 2006. 

[29] K. R. Rajagopal and L. Tao. Mechanics of mixtures, volume 35 of Series on 
Advances in Mathematics for Applied Sciences. World Scientific Publishing Co. 
Inc., River Edge, NJ, 1995. 

[30] C. Truesdell and W. Noll. The non-linear field theories of mechanics. Springer- 
Verlag, Berlin, third edition, 2004. Edited and with a preface by Stuart S. 
Ant man. 

[31] P. Yue, C. Zhou, and J.J. Feng. Spontaneous shrinkage of drops and mass 
conservation in phase-field simulations. J. Comput. Phys., 223(1): 1-9, 2007. 



34 



